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Abstract 

This report documents the development, validation, and applica- 
tion of a spectrally accurate boundary-layer code, WINGBL2, which 
has been designed specifically for use in stability analyses of swept- 
wing configurations. Currently, we consider only the quasi-three- 
dimensional case of an infinitely long wing of constant cross section. 
The effects of streamwise curvature, streamwise pressure gradient, and 
wall suction and/or blowing are taken into account in the governing 
equations and boundary conditions. The boundary-layer equations 
are formulated both for the attachment-line flow and for the evolving 
boundary layer. The boundary-layer equations are solved by marching 
in the direction perpendicular to the leading edge, for which high-order 
(up to fifth) backward differencing techniques are used. In the wall- 
normal direction, a spectral collocation method, based on Chebyshev 
polynomial approximations, is exploited. Spectral accuracy is advan- 
tageous in that 1) the solution is highly accurate even for relatively 
coarse grids, 2) the boundary-layer profiles and their derivatives are 
extremely smooth, and 3) interpolation to other grids can be accom- 
plished with virtually no loss of accuracy. The accuracy, efficiency, and 
user-friendliness of WINGBL2 make it well suited for applications to 
linear stability theory, parabolized stability equation methodology, di- 
rect numerical simulation, and large-eddy simulation. The method is 
validated against existing schemes for three tests cases, including in- 
compressible swept Hiemenz flow and Mach 2.4 flow over an airfoil 
swept at 70° to the free stream. 



1 Introduction 


In the current global economy, the economic pressures on aircraft manufac- 
turers are such that a single major design error could precipitate the financial 
collapse of an entire company. In such a competitive environment, the need 
for accurate, dependable, user-friendly, and robust numerical tools for anal- 
ysis and design cannot be overemphasized. 

In particular, the greater emphasis on drag reduction, transition predic- 
tion, and laminar-flow control technology highlights a current need for highly 
accurate boundary- layer codes. Although it is possible and might seem de- 
sirable to solve the full Navier-Stokes (NS) equations to obtain boundary- 
layer profiles for stability analyses, in practice, NS computations become 
prohibitively expensive whenever the boundary layer is highly resolved. (The 
reader should be aware that stability analyses require not only accurate de- 
termination of the mean flow but also of the first and second derivatives of 
mean-flow quantities; hence, accuracy requirements are stringent.) Conse- 
quently, current state-of-the-art techniques rely on boundary-layer calcula- 
tions that are coupled to outer (inviscid) solutions obtained from the Euler 
equations. For the inner (boundary-layer) solution, spectral methods are 
particularly well suited. In practice, spectral accuracy guarantees extremely 
smooth boundary-layer profiles and ensures thal derivatives and interpolants 
can be computed with no significant loss of accuracy. 

This report documents a spectrally accurate boundary-layer code called 
WINGBL2, which was developed especially for application to analyses of 
crossflow instability and transition on highly swept wings (the dominant 
instability modes if the sweep angle is sufficiently large). In particular, 
WINGBL2 provides the boundary-layer profiles (and their derivatives) that 
are needed for linear stability analyses, for analyses based on parabolized 
stability equation (PSE) methodology, and for the base flow for direct nu- 
merical simulations (DNS) or large-eddy simulations (LES). Currently, we 
consider a quasi-three-dimensional geometry it which the wing is assumed 
to be infinitely long in span. The effects of surface curvature, streamwise 
pressure gradient, and wall suction and/or blowing are incorporated into the 
governing equations and boundary conditions. Because the code is intended 
for use by PSE, DNS, and LES, which consider mean-flow nonparallelism, 
particular attention is paid to the accurate determination of the wall-normal 
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velocity. Procedures for interpolation and extrapolation to other grids are 
also carefully addressed. The code is equally applicable to incompressible 
or compressible flows; however, because of the unnecessary participation of 
the energy equation in the incompressible limit, some performance penalty 
is incurred. 

In the next section, we present the coordinate system and the governing 
equations for the infinite swept-wing geometry. In section 3, we introduce the 
boundary-layer transformations that lead to the transformed governing equa- 
tions. The attachment-line boundary layer is considered as a special case. 
The numerical method is discussed in section 4. Section 5 treats evaluation 
of the wall-normal velocity. Data input and output are discussed in section 
6 , as are data interpolation and extrapolation to other grids. In section 7, 
we close with several validation test cases, among them incompressible swept 
Hiemenz flow and compressible flow over a highly swept wing. 


2 Coordinate System and Boundary-Layer 
Equations 


Let (x*, y*, z*) be a body-oriented orthogonal coordinate system on a swept 
wing, where x * is the arc length from the attachment line in the direction 
perpendicular to the leading edge, y* is the wall-normal coordinate, and z* 
is the spanwise coordinate parallel to the leading edge, as shown in Figure 1 . 
(Throughout this report, dimensional quantities are denoted by an asterisk.) 
The body-fitted system is imbedded in a rectangular Cartesian coordinate 
system, which we denote by (/*, r* , z*) and which will later be used to define 
the wing geometry. Consider the flow of air over the surface of the wing, 
where p*, p*, T*, p*, and k* denote the pressure, density, temperature, vis- 
cosity, and thermal conductivity of the air, respectively. Furthermore, let 
[u*,u*,u?*] T denote the velocity vector in the (x*, j/*, z*) system. 

The three-dimensional compressible boundary-layer equations are given 
in dimensional units for body-fitted coordinates on page 220 of reference [ 1 ], 
Here, we impose two additional assumptions: the flow is laminar, and the 
span is infinite. The latter assumption implies that all terms that contain 
spanwise derivatives (^ 7 ) vanish. The following governing equations result: 
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Continuity: 


ir- {p ’ u "> + a ? (? ' ,V) = 0 


(i) 


^-momentum: 


P^_ du m du* = _ldf_ d_ / ,du*\ 

q dx* P V dy* q dx * + dy * \ dy* ) 

^-momentum: 

ftSdw* , m dw* d ( ,dw*\ 

9 dx * P V dy * dy* dy* / 


( 2 ) 

(3) 


Energy: 


p*u* dT* , ,dT* 
<7 dx* P V dy* 



(4) 


State: 

f = (5) 

where Pr = is the Prandtl number, and fi* s the (ideal) gas constant. 

In the equations above, one metric quantity, q = 1 - y*j^, arises from 
streamwise curvature, where <f> = tan" 1 is the angle in the" 7 ( x *, z*) plane 
of the surface tangent to the body. In boundary-layer theory, 


in which case we obtain 



( 6 ) 


dx* dx* PeU *dx* 


(7) 


for the inviscid flow. Throughout this paper, the subscript e denotes boundary- 
layer edge values. 
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3 Transformed Boundary-Layer Equations 


Transformation of the boundary-layer equations is desirable for two reasons. 
First, Eqs. (l)-(4) are singular at x* = 0. Second, a clever transformation 
will remove most of the timelike evolution in the x* direction, which facil- 
itates marching solutions. In accordance with reference [1], we adopt the 
transformation 


Z = x 


v = 




* 



(8) 

( 9 ) 


where L* is the boundary-layer length scale. defined as 


L* = 



( 10 ) 


and 



(ii) 


is the equation of state under the boundary-layer approximation (Eq. (6)). 

Three elements of the Jacobian of the transformation are readily accessi- 
ble from Eqs. (8)-(9); these are 


and 


dc 

dx* 

d£_ 

dy* 

d r) 
dy* 


= 1 
= 0 
1 

: L*9 


( 12 ) 

( 13 ) 

(14) 


Because the product of the Jabobian matrices of the forward and inverse 
transformations must be the identity matrix, we also ascertain 


dx* 

d£* 


= 1 


(15) 
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and 

dy drj dy* 

dx* dy* d(* 

From Eq. (16), we obtain 

(17) 

y (C,v) — L*(C) f 0(c,r])dT] 

J o 

(18) 

and 


q{t\r]) = 1 - L\C)~ f Q V 0(C,r])dr] 

(19) 

3.1 Transformed Equations for u* > 0 


Provided that neither u* e nor w* vanish, we can define dimensionless velocities 
as lollows: 

F= — 

K 

(20) 

G = — 
w* 

(21) 

V* 

v = 

K 

and 

where Re L is the Reynolds number based on T*; that is, 

(22) 

(23) 

Re, = 

(24) 


Note that Re L - y/W x . The following transformed governing equations 
result from Eqs. (l)-(4), transformation (8)-(9), and definitions (20)-(23): 

Continuity: 


r dF dv F 



x-momentum: 


q 8 (* drjq di) 


z-raomentum: 


Energy: 

CF 86 „86 8 

\-V 

q d£* dr] dr) 

where 


CF8G t 8G 8 
+ V—~ — 


q d(* 8 r) 8 t) [ dr) J 


_8G 


,8F 


= 0 


= 0 


(26) 


(27) 


// 86 
Pr 8q 


-(7-1 )M, 


-2 - (WY, _ 


wi,A d -£\ =o 


8r) 


„ !•** 

0,0 = 1 dt 

Cdni c d( P ;u;) 

° 6* d(* P>*e dt* 

C d K 


fh = 


K d(* 


M, 


u„ 


ii,e 


M - — 

1V1 W,t — 

a? 


\]iRg*T; 


and 


A = ^ m ^ = ^*CC*) 

Me 


(28) 

(29) 

(30) 

(31) 

(32) 

(33) 

(34) 


and where 7 is the ratio of specific heats. For this work, we model viscosity 
(and thermal conductivity) by Sutherland’s law, namely, 


0 3/2 (l + C) 198.6R 

" w= 0 + c ' c = ~tT 

Other viscosity laws are readily implemented. 


(35) 


Although all three velocity components appear above, the governing equa- 
tions remain two dimensional in the sense that the flow variables are functions 
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(£ >7) on ly- The reader should be aware that the transformed continuity 
equation presented on page 397 of reference [1] is in error. 

If the sweep angle xj> = 0, then w* vanishes; however, Eqs. (25)-(28) 
remain valid if G = 0. 


3.2 Attachment-Line Equations 

When the body is blunt, the flow stagnates along the attachment line £* = 0. 
By definition, u e = u =0 for £* = 0, in which case the i-momentum equa- 
tion vanishes. To borrow from the Frobenius method for singular differential 
equations, we differentiate Eq. (2) with respect to x*\ by then redefining F , 
V, and L* as 


F -ll. r- du '- r d < 

J dx *’ Je ~dx * 


(36) 

v - WgV 
f:L*e 


(37) 

and 



r* 1 Me 

“Vrt/r 


(38) 

we obtain an equation system formally identical to that of Eqs. (25)— (28), 
with the exception that the timelike derivative terms vanish and certain 
coefficients differ. Thus, the governing equations for £* = 0 and for £* > 
0 can be put in the following universal form, valid for all £* with proper 
interpretation of F, V, and L*: 

Continuity: 

CdF dV C cl F C c2 

... + a + — -ve = o 

q df Oq q q 


(39) 

r-momentum: 



VFSF dF C A , a r.8Fl 

« ae + K d, "T ( ,_ ^rwj =0 


(40) 

^-momentum: 

CFdG dG a L 

q d(* + drj dq Y dq j “° 


(41) 
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Energy: 


CF dO v <M_d_ 
q d£* + dr) dt] 


fi dd 
Pr dr) 


Cel A 


(sfV t 

(07) “ c ' 2 "l 


. !9G' 

dr) j 


If the body has a sharp leading edge or if £* > 0, then 


and 


and 


= 0 (42) 


Cel — + A)) 

(43) 

o 

$ 

ii 

CM 

o 

(44) 

C„ = A 

(45) 

c.x = (7 - ml. 

(46) 

c. 2 = (7 - mi. 

(47) 

int and £* = 0), then 


Cd = 1 

(48) 

C c 2 - CHq 

(49) 

Cxi = 1 

(50) 

Cel =0 

(51) 

C e2 = ( 7 -l)iV< e 

(52) 


3.3 Boundary Conditions 

The appropriate wall boundary conditions for velocities are 

F(0) = G(0) = 0; V (0) = V wa jj (53) 

where V^ a j] is negative for wall suction or positive for blowing. For the 
temperature at the wall, we consider either the isothermal condition 

0(0) = ^ (54) 
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or the adiabatic-wall condition 



At the far-field boundary 77 = 7/niax 


( 55 ) 


F = G = e = 1 (56) 

In the next section, a transformation will be introduced to eliminate the 
necessity of a far-field condition for V. 


4 Numerical Methodology 


Details of the spectral collocation method can be found in Pruett and Streett 
[2] and Pruett [3]. Here, we provide only a sketch of the method. 

ihe governing equations (39)-(42) are solved by marching downstream 
in the timelike variable Timelike derivatives of the form f^, for example, 
are approximated by high-order backward finite differences. The present code 
allows for mth-order differencing, where 1 < m < 5. In practice, boundary- 
layer solutions are relatively insensitive to the order of the marching scheme; 
we recommend m = 2 or m = 3. The preseni scheme also permits equal 
or unequal steps in . For the initial condition” (£* = 0), the timelike 
derivative terms vanish. For subsequent steps, the order of the scheme is 
limited initially by the number of previous stations available. That is, for 
step 1, m = 1; for step 2, m = 2, and so on, until the full order m of the 
scheme is established. 

At each marching station (r = 0, 1,..., Af), the solution procedure is 
equivalent to solving a system of four coupled nonlinear ordinary boundary- 
value problems in the variable tj. As for finite-difference methods, we parti- 
tion the domain such that 0 = rjo < rji < ... < t]n = »?max> and we define 
the vectors 6 , F , G, and V , whose elements are the discrete approximations 
to their continuous counterparts at the (N + 1) grid points. Moreover, we 
define the discrete differentiation operator [D], whereby, for example, [D]$ 
approximates the continuous derivative Spectral accuracy is achieved by 


10 



approximating the continuous variables with finite expansions in the Cheby- 
shev polynomials 

Tj(fj) = cos [j cos -1 fj (57) 

which are orthogonal on [- 1 , + 1 ] . The natural (Gauss-Lobatto) set of collo- 
cation points associated with Eq. (57) is 

*?i = cos^ 0" = 0, 1, 2, N) (58) 

for which the interpolating polynomial is uniformly bounded as N — ► oo. 
The differentiation operator [£)] is obtained by explicitly differentiating the 
Chebyshev polynomials at the collocation points. Unlike in standard finite- 
difference techniques, [Z)j is a dense, rather than banded, ( N + 1) x (TV + 1) 
matrix. The advantage of spectral methods is that as N — + oo the truncation 
error diminishes faster than N~ p for any finite power P, a property known as 
spectral convergence. An additional advantage of spectral methods is that, 
unlike finite-difference methods, boundary points require no special treat- 
ment. Currently, we compute discrete derivatives by matrix multiplication, 
although a considerable performance gain can be realized for large values of 
N by fast Fourier transform methods. 

For collocation methods, the discrete governing equations are required to 
be exactly satisfied at the collocation points. The discrete governing equa- 
tions are analogous to their continuous counterparts (Eqs. (39)— (42)) ; for 
example, terms like become [ D]V . For initial guesses 0 °, F°, and so on, 
the four discrete governing equations define a residual vector r° of length 
4(iV + 1). The residual iterates r*, where k denotes the iteration index, are 
driven toward zero by Newton’s method. The Jacobian for Newton’s method 
is a nearly dense [4 (N + 1)] x [4(N +1)] matrix. For i > 0, the starting val- 
ues for the iteration are obtained from the converged solution of the previous 
marching step. For spectral methods based on Chebyshev polynomials, Jaco- 
bians are typically quite ill-conditioned. Typically six to seven iterations are 
required to drive the norm of the residual vector to 10 -7 , a level at which the 
solution is smooth to nearly machine precision. To obtain convergence of the 
iteration, all dependencies must be represented by the Jacobian, including, 
for example, the variation of q with 9 implicit in Eq. (19). 

In practice, before the discrete equations are solved, we make one further 
transformation to Eq. (39). By multiplying the continuity equation by q and 
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integrating both sides from r } J _ 1 to Tjj, we obtain 
ov l„ - qV + r 

J Vj~ l 

Integration of Eq. (39) has two advantages. First, after integration, no 
boundary condition is needed for V N (which is now part of the solution). 
Second, a term produced in the integration by parts of exactly cancels the 
term that contains C c2 in Eq. (39). The discrete integration that corresponds 
to the integral in Eq. (59) is also performed by the method described in [ 3 ], 
which makes use of a spectrally accurate quadrature operator [01 in the form 
of an TV x (TV + 1) matrix. 

One final practical consideration requires the use of a continuous transfor- 
mation 77 ( 77 ) from the computational domain [- 1 ,+ 1 ] to the physical domain 
[0,??max]- A well-designed transformation also serves to concentrate points 
near the wall, as is desirable for boundary-layer calculations. The reader is 
referred to reference [3] for details of the specific transformation. In practice, 
the metric dfj/dr), computed either analytically cr numerically to spectral ac- 
curacy, is imbedded directly into the differentiation and quadrature operators 
[D\ and [Q], respectively. 


dF 1 

c ci F + r~ 
dt 


d*l = 0 (j = l,...,N) (59) 


5 Spectrally Accurate Extraction of Wall- 
Normal Velocity 


Historically, the wall-normal v velocity has been the stepchild of boundary- 
layer approximations. Until recently, most stability analyses for boundary- 
layer flows assumed the mean flow to be locally parallel, in which case 
the wall-normal velocity was disregarded altogether. With the advent of 
multiple-scales analyses, PSE methodology, and spatial DNS, each of which 
considers mean-flow nonparallelism, the wall-normal velocity has assumed 
greater importance. This is particularly true for high-speed flows, in which 
nonparallel effects can be quite important. When desired, the wall-normal 
velocity is usually extracted from the continuity equation in physical space; 
this process requires interpolation, which can result in poor accuracy. In ac- 
cordance with Pruett [3], a better approach is to extract v directly from the 
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transformations (23) and (37), in which case the continuity equation remains 
as a check. 

From Eqs. (23) and (37), respectively, we obtain 


v = 


Rei 

w* 


v - 


CF dq 


(C 1 0 ) 


q dx*\ 

VO (£* = 0, blunt body) 


( 60 ) 

( 61 ) 


The extraction of v is straightforward for the attachment-line problem 
£* = 0. However, for £* ^ 0, Eq. (60) requires an expression for J^-. By 
explicitly differentiating transformation (9) with respect to x*, we obtain 


0 


dr] 

dx* 



J/l I 
2 L* [t* 


1 d{Rei) 

Rei d£* 


(62) 


where Re i = p* e u* e /]i* e is the edge unit Reynolds number. To obtain Eq. (62), 
we have used Leibniz’ Rule to differentiate Eq. (9) beneath the integral sign. 
The integrations necessitated by the first term in Eq. (62) and y * (see Eq. 
(18)) in the second term of Eq. (62) are computed to spectral accuracy by 
the integration procedure described in detail in reference [3]. 


6 Code Input and Output 

6.1 Input 

Four input files are required to define the physical parameter values, the pa- 
rameters of the method, the geometry, and the boundary-layer edge data. 
Physical constants are set in a file “constants. h,” which is incorporated into 
subroutine “econl.f” by a FORTRAN include statement. This file normally 
resides with the source code. Default values are supplied for any physical con- 
stants not specified. Parameters of the method are read from standard input 
(FORTRAN unit 5) in the initialization subroutine “init.f.” We have named 
the standard input file “wingbl2.d.” The wing geometry is specified in file 
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“geom.dat” and read by subroutine “sgeoml.f,” which uses cubic splines to 
obtain accurate interpolated values for f, <j > , and so on. The boundary-layer 
edge data are specified in file “edge.dat” and read by subroutine “econl.f,” in 
which splines are also exploited for accurate interpolation of the edge data. 
Sample listings of the input data files are provided in appendix A. These 
files and the source code are documented such that the role of most param- 
eters is self-explanatory. In the following subsections, we provide additional 
documentation for those parameters whose usage is not self-evident. 


6.2 Determination of Arc Length 

File “geom.dat” contains the chord length c* (CHORD), the sweep angle xb in 
egrees (SWEEPD), a parameter KSPLIT, whose function will soon become 
evident, and tabulated values of l*/c* and r*/c‘, which define the wing ge- 
ometry relative to some origin whose precise coordinates are unimportant. If 
the parameter CHORD is negative, the negative sign is interpreted to mean 
that the^ chord length c* (normal to the leading edge) has been specified in 
ieu of c (chord length parallel to the free-stream velocity vector). On the 
basis of Figure 1, the relationship between these two values is c* n = c * cos xb. 

For the case (CHORD < 0), the tabulated values l* and r* are assumed to 
be normalized by c*. 

Considerable care was devoted to ensure that geometric and edge data 
are evaluated accurately. Specifically, to avoid singularities in the calculation 
of arc length, we employ the two equivalent expressions 



where Eq. (59) is used in the attachment-line region and Eq. (60) is used 
elsewhere. The switch from Eq. (59) to (60) is made at grid point KSPLIT 
(an input parameter referred to above), typically near the point at which 
dr* /dr - 1 . The integrands on the different regions are independently fit 
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with splines, and the resulting spline functions are integrated analytically. 
To ensure continuity of first and second derivatives at the interface, the 
end conditions on the matching splines are adjusted as follows. It can be 
determined from elementary calculus that 


dr * 1 

dF = dFJdF 


(65) 


and 

dV / dr*\ 3 (Pi* 

dl * 2 “ “ { dl* j dr* 2 ( 66 ) 

An initial guess is made for as an end condition for one spline; the 
other end condition is determined from Eq. (66). The assumed value is 
subsequently adjusted with the secant method until Eq. (65) is satisfied. At 
convergence, functional values and first and second derivatives are consistent 
at the interface. 


Similarly, <f) and are obtained from two different sets of expressions: 

* = ATAN2^,1.); ^ = C os 3 ^ (67) 

^ ATAN2 < 1 "^ < 68 > 

where ATAN2 is a particular two-argument FORTRAN arctan function. 


6.3 Gas Relations for Edge Conditions 


Our interest is in high-efficiency wings without shocks. Consequently, we 
assume that the flow is isentropic and that the fluid is a perfect gas. (It 
is straightforward, however, to modify the gas relations to consider non- 
isentropic flow.) From file “edge.dat,” we obtain the free-stream Mach num- 
ber Moo (MINF), the pressure p (PINF), and the temperature (TINF), 
as well as tabulated values of the chordwise coordinate /*/c* (XC) (scaled 
either by c* or c*), the pressure coefficient Cp e (PC), and the suction coeffi- 
cient Cq e (QC). From this data, the total free-stream velocity (UTINF) 
is computed as 

VI = M„yJ-,Rg-T^ (69) 
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Edge pressure data are extracted from the pressure coefficients (PC) as 


Cp e = 



From Eq. (70), we obtain 


~r = i + \c Pc Ml 
2 


( 70 ) 


(71) 


From the isentropic perfect gas relations given in reference [4], we obtain 



(72) 


and 


U* 2 2 

(—) = 1 + 

+ (7-l )Ml 



(73) 


For the infinite swept wing, the spanwise edge velocity is constant, namely, 


W *e=UZo S ™i> (74) 

From Eqs. (73)-(74) and the expression 

U* = \f< + K 2 (75) 

we obtain u*. The contribution of the wall-normal velocity to the total 
velocity is assumed to be inconsequential. Given these quantities, the vari- 
ations along the boundary-layer edge of all other quantities of interest (e.g., 
Reynolds number) readily follow. Two additional parameters head file “edge.dat.” 
These are IFIRST, which relates the first indices of the tabulated values of 
files “geom.dat” and “edge.dat,” and ISTAG, which identifies the index of 
the stagnation point (attachment line). At present, the tabulated locations 
l*/c* of “edge.dat” are assumed to be a contiguous subset of the values listed 
in “geom.dat.” Furthermore, the stagnation point is assumed to be known. 
However, the spline interpolation algorithms are quite general, and these 
restrictions are easily relaxed. 
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6.4 Wall Suction 


From the tabulated suction coefficients Cq e (QC), we derive values for 
as follows. From the definition 


Cq e = 



( 76 ) 


we obtain either 


for u* / 0 or 


I wall — 




Cq. 


q I wall “ 



(77) 

(78) 


for attachment-line flow. From Eqs. (77) and (78) above and Eqs. (23) and 
(37), appropriate expressions for are readily obtained. (Note that Eq. 

(23) simplifies considerably at the wall because vanishes there.) 


6.5 Output 


The output consists of a file “legend” and a collection of numbered files begin- 
ning with “sta_0000” and numbered with increments of IOUT (an output pa- 
rameter), all placed in subdirectory “sta” (for “stations”). Each “sta_xxxx” 
file corresponds to a particular value of x* and contains the dimensionless 
quantities 6, F , G, and v and their first and second derivatives with respect 
to y. Appended at the end of each “sta_xxxx” file is a list of approximately 
30 relevant dimensional and dimensionless values that apply at the given 
station, among them x *, local Mach number, Reynolds number, displace- 
ment thickness, and boundary-layer edge conditions. The file “legend” tab- 
ulates summary information from the “sta_xxxx” files. Sample “legend” and 
“sta_xxxx” files are given in appendix B. Lengths in the output files are non- 
dimensionalized at the user’s option (ISCALE) either by the boundary-layer 
length scale L* (Eq. (10)) or by the boundary-layer displacement thickness 
8*, which is defined as 


- 1 ° 

Jo 


1 - 




PeK 


/imai 

dy* = L* I [6 — F] df] 

Jo 


(79) 
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Among the quantities tabulated in the “legend” file are 8* and the ratio 8* / L* 
as functions of x* . As the Mach number increases, the ratio 8* J L* increases 
dramatically from its value of about 1.7 at M ^ = 0. Consequently, 8* is the 
more practical length scale, although L* is the conventional scale in stability 
theory. A recommended rule of thumb is to adjust r) max (ETAJV1AX, an 
input parameter) to ensure that the computational domain is at least 3£* in 
thickness. Otherwise, the outer boundary condition will “pinch” the solution. 

In addition to the output described above, several diagnostic files are 
produced; among these is the standard output file to FORTRAN unit 6. 


6.6 Interpolation and Extrapolation 


Because the present code is designed specifically for applications to stability 
analyses, DNS, and LES, we have paid special attention to the process of 
interpolating and extrapolating data to other grids. For this purpose, a 
companion code INTBL has been written. In brief, the procedure is as 
follows. In the output from WINGBL2, lengths are scaled either by L* or 
8\ each of which grows with x*. The outer edge of the boundary layer 
typically lies between 2 and 4 displacement thicknesses from the wall. Let 
ye{% ) denote the location of the boundary-lay ir edge, where y = y* J L* 
To interpolate to a grid that is uniform in y , for example, we first perform 
spectrally accurate interpolation in the interior region (0 < y < y e ) coupled 
with analytic extrapolation outside the boundary layer (y > y e ). This step 
is followed by high-order (typically fifth) polynomial interpolation in x *. 

Analytical extrapolation of u*, T * , and w* is trivial because, for example, 
u = u e i n ^he outer flow. Extrapolation of w*, however, requires some care. 
From the continuity equation (Eq. (1)), we derive the following ordinary 
differential equation for the far-field behavior of v: 


(1 - a oy)j7 ~ a ov = -/? 

where a 0 and v are defined in Eqs. (29) and (22) , respectively, and 

L* d(/>X) 

p*u* dx m 


(80) 


(81) 
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Equation (80) has the exact solution 


v(y) = v e 


1 - ftp ye 
1 - <*oJ/ 




1 -a 0 y 


where v e = v(y e ). From Eq. (82), we also obtain v'(y ) and v"(y). 


(82) 


7 Code Validation 


We present three test cases that collectively validate WINGBL2 for compress- 
ible and incompressible flows, for two-dimensional and quasi-three-dimensional 
flows, and for flows over bodies with and without streamwise curvature and 
with and without streamwise pressure gradients. 


7.1 Two-Dimensional Compressible Flow Along a Flat 
Plate 

Here, we consider a two-dimensional Mach 4.5 flow of air along a flat plate 
with a sharp leading edge. For a thin plate and a thin boundary layer, the 
shock is oblique and weak and the flow can be considered to be approximately 
isentropic. The parameters of the flow are = 4.5, T , = 110°R, Re x = 
2.888 x 10 6 , 7 = 1.4, and Pr = 0.7. The pressure and the suction coefficients 
are zero. The geometric parameters are xp = 0, with a sharp leading edge 
(IBLUNT=0) and no streamwise curvature (r*(/*) = 0). In the absence of a 
streamwise pressure gradient, a similarity solution exists. Figure 2 compares 
results from the present code with those obtained from the spectrally accurate 
boundary-layer code of Pruett and Streett [2]. Results from the two methods 
are indistinguishable. 


7.2 Incompressible Swept Hiemenz Flow 

The incompressible swept Hiemenz problem, which was studied thoroughly 
in Malik et al. [5], represents an idealized stagnation-region flow on a swept 
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flat plate. In the absence of body curvature, u * varies linearly with x* (i.e., 
du*/dx* is constant), and a similarity solution exists for which the boundary- 
layer thickness remains fixed with x*. The reader is referred to reference [5] 
for details of the problem and the geometry. Figure 3 presents the results 
obtained by the present method with those of reference [5] at a station x* 
for which Re^ — 500 and u* e = w* e . The results of the two methods are 
indistinguishable. Results for incompressible flow were obtained from the 
present code by setting = 0.01, p*^ = 2000 psf, and = 520°R, and 
by adjusting the velocity gradient du*/dx* and the sweep angle ip to obtain 
Rcl = 500 and u* = w* at some arbitrary value of x*. The parameters 
selected were xp = 50.7° and du*/dx* = 0.0961 sec -1 , which gave x* = 90 
ft as the station of interest. We point out that for the present method 
the algorithm actually switches from the attachment-line equations to the 
equations of the evolving boundary-layer flow; however, the solutions are 
virtually identical, as they should be for this unique problem. 


7.3 Flow on a Supersonic, Highly Swept Wing 


Here, we compute the boundary layer on the upper surface of an infinitely 
long shockless wing with a chord length 20 ft; the wing is swept 70.22° to 
a free-stream velocity of Mach 2.4. The input files presented in appendix 
A contain the data for this test case (condensed for brevity). This case 
represents a stringent test for boundary-layer codes. Streamwise curvature 
in the attachment-line region is quite strong, and regions of both favorable 
and adverse pressure gradients exist. The maximum Mach number attained 
in the direction perpendicular to the wing is cbout 0.95, in which case a 
shock is avoided and isentropic gas relations apply. Figures 4 (a,b,c) compare 
the present results with those obtained from the code of Wie [6] at stations 
x * = 6.9xl0 -3 , 0.851, and 4.027 ft, respectively. For both the finite-difference 
and spectral methods, approximately 40 grid points were used across the 
boundary layer. The agreement is good at all three stations, except for some 
disagreement in the wall-normal velocity. Both solutions fail to converge for 
l*/c* > 0.6, where there is a strong adverse pressure gradient. With strong 
suction in the trailing-edge region, the present code predicts that attached 
flow can be maintained. To obtain a solution for l*/c* > 0.6 without suction, 
interaction must be allowed between the invisc d region and the boundary 
layer. In Figure 5, we compare second derivatives of the u - velocity component 
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at station x* = 4.027 ft. The agreement is quite good. 

For applications to stability analyses, we prefer to use approximately 100 
points across the boundary layer, even for the spectral method. Such fine 
resolution is unnecessary but desirable for two reasons. First, for the present 
spectral method with N = 100, the velocity and temperature profiles and 
their first and second derivatives are smooth to 13, 10, and 7 significant digits, 
respectively. For stability analyses, the smoothness of the profiles and their 
derivatives, in addition to accuracy, is important. Second, a high-density 
grid ensures virtually no loss of accuracy in interpolation to other grids. 

In summary, the present method is well suited for stability analyses that 
use any of the following numerical tools: LST, PSE, DNS, or LES. The spec- 
tral accuracy of the present scheme and its accurate determination of the 
wall-normal velocity render WINGBL2 particularly attractive for applica- 
tion to high-speed boundary-layer flows, which are quite sensitive to subtle 
changes in the mean state. 
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Appendix A: Input Files 


Input file “constants. h” 


rgas = 1716. 
col = 2.27e-8 
co2 = 198 . 6 
pr = 0.72 
gamma = 1.4 


Input file “wingbl2.d” 


1 / I BLUNT = 0 (sharp leading edge); 1 (blunt body) 

0 /ISCALE = 0 (scale lengths by L*); 1 (scale by delta*) 

60 /NSTEPS = number of steps in marching (xi) direction 
-0.05 /DX = dimensional increment in x if > 0 (otherwise use input grid) 

42 /NP = N+l = number of grid points in wall-normal (eta) direction 

0 /ISTART = start switch: 0 (fresh start); i > 0 (restart from i) 

3 /M = max. order of backward difference method [2 or 3 RECOMMENDED] 

1 /IOUT = output stride for sta_xxxx files 

5.e-7 /RLIM = maximum residual allowed [KEEP DEFAULT 5.e-7] 

20 / ITMAX = iteration limit [KEEP DEFAULT 20] 

0.7 /STR = parameter of tanh stretching [KEEP DEFAULT 0.7] 

12.0 / ETA_MAX = outer bound in wall-normal variable eta 

-.35 /RO = initial guess for F = 1 - exp(ro*eta) [KEEP DEFAULT -.35] 
300. /TWALL = dimensional wall temp. (Rankine--irrelevant for iadiab=l) 
1 / IADIAB = wall boundary condition: 0 (isothermal); 1 (adiabatic) 

notes: 1) to run, create the subdirectory /sta 

2) input files needed 

'geom.dat' geometry data 
'edge.dat ' boundary- layer edge data 

'wingbl2.d' standard input containing parameters of 

method 

3) for cubic spline interpolation, edge data file must have 
at least 4 points 

4) at completion, subdirectory /sta will contain 
file 'legend' 

nsteps/iout numbered files / sta_xxxx / 
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Input file “geom.dat” 


20 /number of descriptive header records 
# 

####### ########################################################## 

# 

# title: SUPERSONIC SWEPT-WING CASE 

# 

# airfoil: NACA16006 

# 

# geometry data: 

# 


# 

# 

# 

# 

# 

# 


CHORD: chord length (c*) if > 0 

or normal chord (c_n*) if < 0 
SWEEPD: sweep angle (psi) in degrees 

KSPLIT: index at which to switch from 

Eq. (63) to Eq. (64) 


tabular values of length 1* and r* normalized 

by c* or c_n* according to value of CHORD 


################################ ################################# 

# 

20.0000 /chord 

70.2180970318 /sweepd 
8 /ksplit 

2 . 99717 07E-3 -1 . 0969796405742E-3 

1 . 8810658E-3 -8 . 7168943307489E-4 

9 . 4530382E-4 -6 . 2043138039564E-4 

2 . 6244059E-4 -3 . 2962356850761E-4 

0 . 0 . 

2 . 6244059E-4 3 . 2962356850761E-4 

9 . 453 0382E-4 6 . 2043138039564E-4 

1 . 8810658E-3 8 .7168943307489E-4 

2.99717 07E-3 1 . 0969796405742E-3 


0.9824749 1 . 0061736003393E-3 

0.9865789 7 . 8976675042668E-4 

0.9903595 5 . 7946564411269E-4 

0.9938442 3 . 7687616857631E-4 

0. 9970505 1 . 8348090931462E-4 

1 . 0 . 
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Input file “edge.dat” 


9 /number of descriptive header records 
# 

########## ####################################################### 

# 

# Freestream data, pressure and suction coefficients 

# 

# tabular values of normalized length 1*, cp_e, cq_e 

# 

###################################################### ########### 

# 

5 /ifirst (index in 'geom.dat' corresponding to first point) 

1 /istag (index of stagnation point/attachment line) 

2.4 /rminf (freestream Mach number) 

149.76 /pinf (freestream pressure) 

390.0 /tinf (freestream temperature) 


0 . 

0.1347191240328 

0 

2 . 6244059E-4 

0.1132312733923 

0 

9 . 453 03 82E-4 

7. 01245573943 67E-2 

0 

1.8810658E-3 

4. 6016691556783 E- 2 

0 

2 . 997 17 07E-3 

3 . 1253217343323E-2 

0 

4 . 27 05112E-3 

2 . 1114260233283E-2 

0 

5 . 6953 584E-3 

1.3796344532341E-2 

0 

7 .2747553E-3 

8.3160311603793E-3 

0 

9 . 0147518E-3 

4.0188934669989E-3 

0 

1 .0925392E-2 

5. 434347772872 IE-4 

0 

1 .301854E-2 

-2 . 14890547270 16E-3 

0 

1.5307413E-2 

-4.103024062921E-3 

0 

1 . 780377 5E-2 

-6.3 03813862427 6E-3 

0 

2 . 0 52 6798E-2 

-8.4553393854651E-3 

0 

2 . 34928E-2 

-9.8847372967642E-3 

0 

2 . 6720202E-2 

-1 . 1098457423 478E-2 

0 


0.9824749 

2 . 358551808154 4E -2 

0 

0.9865789 

2 . 8966689722272E-2 

0 

0.9903595 

3 .5253 108067807 E-2 

0 

0.9938442 

4 . 3951539866815E-2 

0 

0.9970505 

5.3919871571419E-2 

0 

1 . 

6. 594498775923 5E-2 

0 
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Appendix B: Output Files 


Out file file “legend” 


step x* del*/L* delta* Mach_u rmax imax 

0 0 .0000000E+00 0 . 1430674E+01 0 . 2222551E-03 0 . OOOOOOOE+OO 0 . 1350571E-07 6 

2 0 . 1435503E-01 0 . 1763038E+01 0 . 3 173347E-03 0 . 5206726E+00 0 . 2154925E-06 11 

4 0 . 3 125234E-01 0 . 2366497E+01 0 . 5700238E-03 0 . 6863386E+00 0 . 2066538E-06 13 

6 0 . 5126869E-01 0 . 2689383E+01 0 . 8060416E-03 0 .7568990E+00 0 . 3413945E-06 12 

8 0.7498714E-01 0 . 2913725E+01 0 . 1041978E-02 0 .7961337E+00 0 . 44113 12E-06 11 

10 0 . 1030743E+00 0 . 3100542E+01 0 . 1290048E-02 0 . 8208850E+00 0 . 1781216E-06 11 

12 0 . 1362966E+00 0 . 3217949E+01 0 . 1532241E-02 0 . 8375850E+00 0 . 2482159E-06 7 

14 0 . 1755287E+00 0 . 3329082E+01 0 . 1791834E-02 0 . 8520063E+00 0 . 1293438E-06 10 

16 0.2217 67 6E+00 0 . 3434735E+01 0 . 2072884E-02 0 . 8614153E+00 0 . 1185993E-06 9 

18 0 . 27 61402E+00 0 . 3506434E+01 0 . 2356517E-02 0 . 8695486E+00 0 . 1832802E-06 8 

20 0.3399068E+00 0 . 3577373E+01 0 . 2663127E-02 0 . 8760325E+00 0 . 1388041E-06 8 

22 0 .4144507E+00 0 . 3632934E+01 0 . 2982440E-02 0 . 8814904E+00 0 . 2513429E-06 7 

24 0 .5012724E+00 0 .3676531E+01 0 . 3315544E-02 0 . 8863561E+00 0 . 1103803E-06 7 

26 0 . 6019561E+00 0 . 37 18341E+01 0 . 3671124E-02 0 . 8904525E+00 0 . 4261433E-06 6 

28 0 .7181306E+00 0 . 3756340E+01 0 . 4047261E-02 0 . 8942131E+00 0 . 8833897E-07 7 

30 0 . 8514025E+00 0 . 3767767E+01 0 . 4416550E-02 0 . 8979295E+00 0 . 8888528E-07 7 

32 0 . 1003277E+01 0 . 3779752E+01 0 . 4805359E-02 0 . 9018746E+00 0 . 3266885E-06 5 

34 0 . 1175043E+01 0 . 3792276E+01 0 . 5213497E-02 0 . 9055657E+00 0 . 8976225E-07 6 

36 0.13 67 650E+01 0 . 3795478E+01 0 . 5624686E-02 0 . 9094002E+00 0 .2393598E-06 6 

38 0 . 1581563E+01 0 . 378993 1E+01 0 . 6034509E-02 0 . 9135010E+00 0 . 4500929E-06 6 

40 0 . 1816618E+01 0 . 3781879E+01 0 . 6447810E-02 0 . 9178708E+00 0 . 4005004E-06 6 

42 0 .2071904E+01 0 . 3774961E+01 ■ 0 . 6867179E-02 0 .9222863E+00 0 . 4239713E-06 6 

44 0 .2345658E+01 0 .3768265E+01 0 . 7287375E-02 0 . 9266923E+00 0 . 6726634E-07 7 

46 0 .2635245E+01 0 . 3762494E+01 0 . 7705623E-02 0 . 93 10882E+00 0 . 326723 6E-06 6 

48 0. 2937203 E+ 01 0 . 3764590E+01 0 . 8133370E-02 0 . 9350821E+00 0 . 1027020E-06 5 

50 0 .3247372E+01 0 . 3767901E+01 0 . 8553584E-02 0 . 9387523E+00 0 . 2202283E-06 5 

52 0.3561097E+01 0 . 3765214E+01 0 . 8945024E-02 0 . 9422244E+00 0 . 1433897E-06 7 

54 0 .3873515E+01 0 . 3783265E+01 0 . 93 68727E-02 0 . 9451995E+00 0 . 3759408E-06 7 

56 0 . 4179840E+01 0 . 3878009E+01 0 . 9974963E-02 0 . 9456833E+00 0 . 1965184E-06 7 

58 0 . 4475640E+01 0 . 4109018E+01 0 . 1094322E-01 0 . 9424782E+00 0 . 1737138E-06 7 
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Outfile file “sta_0058” 


3 /number of header records 
0.447564030 67 253 E+ 01 /x* 

0.10943219324209E-01 /delta* 

0.41090179017183E+01 - . 173 3 0 63 1 1150 0 0E+01 /dol, phi(deg) 

42 

0 . 000000000000000E+00 0 . 10 0000 000 0000 00E+01 0 . 0000 000 000 0000 0E+00 

0 .2049257838884E+01 0 . 0000000 0 0 0 00 0E+00 0 . OOOOOOOOOOOOOE+OO 0 . OOOOOOOOOOOOOE+OO 
0. 1775266507739 E- 11 0 . 1499 8643 9 5883E+00 0 . 1650165890827E+00 0 . 187 04552 05679E-13 
-. 4736774938549 E- 01 0 . 423 5754 057480E-01 - . 1875125248799E-09 0 . 7 17 58 622 11209E-03 
0. 142881152 183 827E-01 0 . 100000128910483E+01 -0 . 526462043 1 19655E-04 

0 .2049253001368E+01 0 . 2147348282272E-02 0 . 23 57777379672E-02 0 . 73 0829887269 IE-07 
677 3 10282386 8E -03 0 . 1505918742467E+00 0 . 1650168731993E+00 0 . 102 1837297 079E-04 
-.4743991377321E-01 0 . 42387393 13721E-01 0 . 4042780866942E-0 4 0 . 7127460511339E-03 


0 .145894283 15580 1E+02 0 . 10013 1628995567E+01 0 . 260143956494611E+01 

0 . 1000000000000E+01 0 . 1000000000000E+01 0 . 1000000000000E+01 0 .4387097619112E-02 
158786406245 8E -08 - . 3 29 69 18293 927E-1 1 - .7673861546209E-11 0 . 4498749479165E-05 
0. 345853 1461220E-08 -. 32459923 63787E-10 - . 7 122808454163E-09 - . 8246321650703E-09 
27 /no. of trailer records 

0 .447564030672532E+01 /x* 

0 .657951999999998E+00 /l*/c* 

0. 95104224477017 8E -02 /t*/c* 

-0.1733 063 11149996E+01 /phi (deg) 

0 . 109432 193242091E-01 /displacement thickness* 

0 .410901790171832E+01 /delta*/L* 

0 . 248885825547478E+01 /edge total Mach no. 

0 . 942478246658577E+00 /edge u Mach no. 

0 . 23 03 50823 102943E+01 /edge w Mach no. 

0 . 63 1016966003 098E+06 /edge unit Reynolds no.* 

0 . 282420496727099E+07 /edge Re no. based on x* 

0 . 168053710678193E+04 /edge Re no. based on L* 

0 . 69053 5705626887E+04 /edge Re no. based on delta* 

0 . 894399863413379E+03 /edge u-velocity* 

0 . OOOOOOOOOOOOOOOE+OO /wall v-velocity* 

0 . 218600000000920E+04 /edge w-velocity* 

0 . 130386856151895E+03 /edge pressure* 

0 . 202694083 069417E-03 /edge density* 

0 . 374865488202295E+03 /edge temperature* 

0 . 768196040225634E+03 /wall temperature* 

0 . 287297442032797E-06 /edge viscosity* 

0 . 140000000000000E+01 /gamma 

0 . 719999999999999E+00 /Prandtl number 

0 . 227000000000000E-07 /first Sutherland constant* 

0 . 198600000000000E+03 /second Sutherland constant* 

0 . 171600000000000E+04 /ideal gas constant* 

0 . 173713816664645E-06 /rmax 
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Figure 2: Comparison of present, results (P) with those of Pruett and Streett 
(PS) for two-dimensional Mach 4.5 flow over insulated flat plate. 
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(a) /*/c* = 0.001. 


Figure 4: Comparison of velocities obtained by present method (P) with 
results obtained from Wie’s (W) code for highly swept wing in Mach 2.4 
flow. 
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